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I . Introduction 

This report describes three procedures for estimating the 
position of an object. Each procedure is based on a model that 
describes the position of an object in terms of the position of a 
point on a plane surface. 

The first procedure is developed in Section II. The 
procedure is a method for combining into a composite estimate two 
or more point estimates of the position of an object that each 
have confidence regions that are bounded by an ellipse. In the 
model on which it is based, the coordinates of each estimate are 
the values of random variables that have a bivariate normal 
distribution with a mean vector whose elements are the unknown 
coordinates of the object. A user of the procedure must specify 
the coordinates of each estimate and the major and minor axes of 
the ellipse that bounds its confidence region. 

The second procedure is developed in Section III. The 
procedure is a method for generating a position estimate from 
bearings that are on or from two or more reference points of 
known position. In the model on which it is based, bearing 
errors are independent normally distributed random variables. 
Based on the model, the coordinates of a position estimate that 
is determined by the procedure are values of random variables 
that have a bivariate normal distribution with a mean vector 
whose elements are the unknown coordinates of the object. A user 
of the procedure must specify the standard deviations of the 
bearing errors. An extension of the procedure is also developed 
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in this section. It is a method for generating a position 
estimate from bearings that are on or from two or more reference 
points for which the position of one or more of the reference 
points is uncertain. 

The third procedure is described in Section IV. The 
procedure is a method for generating a position estimate that is 
determined from lines of position. In the model on which the 
development is based, lines of position are straight lines, lines 
of position that are determined by observation are parallel to 
true lines of position and the algebraic distance between a line 
of position based on an observation and a true line of position 
is the value of a normally distributed random variable. A user 
of the procedure must specify the standard deviations of these 
random variables. 

In the models on which the three procedures are based, 
systematic errors represent mean values that are set to zero. 

This suggests that the procedures should be used only were the 
bias in the measurements that determine bearings or lines of 
position is known or is known to be negligible. 
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II. A Composite Position Estimate 

The procedure that is developed in this section is for 
combining position estimates for an object that are generated 
from data from independent sources. It is based on the following 
model: In a rectangular coordinate system, the coordinates of 

each position estimate are values of random variables that are 
determined by an independent bivariate normal distribution with a 
covariance matrix whose elements are known and with a mean vector 
whose elements are the unknown coordinates of the object. The 
above model implies that the natural logarithm of the likelihood 
function L for a set of n estimates can be expressed as 
follows: log L = K - 3 ~ 2£) where K is a 

constant, x^ is an estimate of the mean vector with components 
x^ and y^ , x is the mean vector with components x and y, the 
unknown coordinates of the object, and 2 ^ is the covariance 
matrix with elements , a^- and o-*. The composite maximum 

likelihood estimates x and y of the unknown coordinates x 
and y make log L a maximum which implies the following 
relations : 



3 (log L) 
5x 



x=x 



0 



and 



3 flog L) 

dy 



x=x 



0 



y=y 



y=y 



After evaluating the derivatives, they can be written as: 



A X + B y = D 



B X + C y = E 
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where : 



A = a^, B = , C = Cj , D = 2" (a^x^ + b^yj , 

E = Z;" (b^i. + c^y^) , a^ = / d- , , 

c. = / d, , d, = a ?5 - a,j* 

The solutions to the equations are: 

X = (CD - BE) / (AC - B^) 
y = (AE - BD) / (AC - B^) 

Since the estimates x and y are linear combinations of the 
estimates Xj and y^ , they are determined by a bivariate 
normal distribution that is defined by its mean vector with 
components x and y and its covariance matrix with elements 

a\, and a-. With E the expected value operator: 

X = E((CD - BE) /(AC - B^) } 

y = E((AE - BD)/(AC - B^) } 

a? = E((CD - BE) - E(CD - BE))V(AC - B^)^ 

= (C^(F+I+2L) - 2CB(H+K+M) + B^ (G+J+2N) )/ ( AC - B^)^ 

a? = E((AE - BD) - E(AE - BD))V(AC - B^)^ 

= (A^(G+J+2N) - 2AB(H+K+M) + B^ (F+I+2L) }/ (AC - B^)^ 

aj- = E( [ (CD-BE) -E (CD-BE) ] [ (AE-BD) -E (AE-BD) ] }/ (AC - B^)^ 

= { (AC+B^) (H+K+M) -CB(F+I+2L) -BA(G+J+2N) }/ (AC - B^)^ 
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where : 



F = S? a^a^- , G = 2 ? , H = 2 ? a b , 1 = 2 ? b^a^- , 



. 2^2 



2^2 



J = 2 ? c?a?- . K = 2 ? b c a^- , L = 2 ? aba-- , M = 2 ? fa c +b^)a -- 



N = 2J b^c.a.j- 



Without loss of generality, the origin of the coordinate 
system can be the point determined by the composite estimate, its 
positive y-axis can be north and its positive x-axis east. In 
Section II of this report, it is shown that the axes of the 
ellipse that bounds a minimum area confidence region for a 
composite estimate are coincident with a primed coordinate system 
which is transformed to this coordinate system by a coordinate 
axes rotation equal to 7 where 7 is defined by the relation: 
tan 27 = 2 aj- /(a? - o?) . 

It can be shown that a minimum area confidence region for a 
two dimensional estimator whose components have a bivariate 

normal distribution is bounded by an ellipse whose major axes are 

• 1/2 

determined by 2 kai- and 2 ka-. where k = [-2 log(l - p) ] 

and x' and y' refer to a coordinate system whose axes are 
parallel to the axes of the bounding ellipse. In this system: 



II 

CM <X 


-’I 


COS^ 


7 - 


2 




sin 


7 


cos 7 


+ 


2 • 2 
aj sin 


= 


“1 


sin^ 


7 + 


2 




sin 


7 


cos 7 


+ 


2 2 

0 ^ cos 
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The model that is described above can be used to find a 



composite position estimate from two or more position estimates 
for an object. To show this, first assume that the estimates 
satisfy the conditions that are required by the model. Then, let 
5, be the direction of the major axis of the ellipse that bounds 
the confidence region that is associated with the i th position 
estimate. Next, orient the primed coordinate system defined 
above so that its x'-axis is coincident with the major axis of 
the ellipse. Then, the elements of the covariance matrix of the 
bivariate normal distribution that determines the estimate are 
given by the following relations: 



2 2 *2 
afx = at;, sin 




+ 




cos^ S 


COS^ 

^ ly ^ix 




+ 


atj. 


sin^ S 


cj... = (a-. - 

ixy ' ^ IX 


o%.) 


sin 


8 , cos 



As an example of the use of the procedure, suppose that each 
bounding ellipse were a circle. In this case, the coordinates of 
each position estimate would be determined by a circular normal 
distribution and, for i = 1 to n , the elements of the covariance 
matrices of the circular normal distributions would satisfy the 
relations: = a, and = 0. The composite estimate 

would be: 

X = (x^ / o^)] / [S; (1 / a^) ] 

Y = [S; (y^ / a^)] / [(2J 1 / a^)] 
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And the elements of the covariance matrix of the distribution 
that determine the coordinates of the estimate would be: 

(7? = n / [SJ (1 / , a? = n /[(Sj (1 / a^)f and = 0 . 

Consequently, the distribution would also be circular normal and 
the bounding curves for minimum area confidence regions would be 
circles . 

A program called COMP that is based on the procedure that is 
developed in this section is listed in Appendix 1. COMP can be 
used to generate composite position estimates and corresponding 
confidence regions. As an example of its use, suppose three 
independent position estimates are given in Table 1. 



i 




Yi 


TABLE 1 


MJj 


MI, 


k 


1 


-3.7 


18.1 


59° 


36 


20 


2 


2 


11.8 


• 

00 


105° 


38 


10 


2 


3 


0 


0 


146° 


50 


24 


2 



In Table 1, MJ is the length of the major axis and MI is the 
length of the minor axis of the ellipse bounding a minimum area 
confidence region. If the distance unit is the nautical mile, 
the coordinates of the composite estimates of the coordinates of 
the object are: x = -2.69 nm and y = 12.41 nm. For k = 2, 
the lengths of the major and minor axes of the bounding ellipse 
of the minimum area confidence region are: MJ = 17.33 nm and 

MI = 8.85 nm and the direction of its major axis is 103.77°. 
The results are shown in Figure 1. 
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North 




Figure 1. The position estimates and their associated confidence 
regions for the example on Page 7. The position estimates are at 
the center of their corresponding confidence region. The numbers 
indicate the order of the estimates in the table on Page 7. The 
composite confidence region is in black. 
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III. A Position Estimate Based on Bearings 

The procedure that is developed in this section is for 
determining a position estimate for an object from a set of 
observed bearings from or on two or more reference points of 
known positions. In the development, reference points are 
referred to as stations. The procedure is based on the following 
model: A set of bearings taken on an object from a station are 

the values of independent normally distributed random variables 
with known standard deviations but with means that are equal to 
the unknown true bearings. The model implies that the natural 
logarithm of the likelihood function L for a set of n 
bearings taken on an object can be expressed as follows: 
log L = K - 1/2 / e^ where K is a constant, 

is an observed bearing, is the unknown true bearing and e^ 

is the known standard deviation. For i = 1 to n , the maximum 
likelihood estimates are the solutions of the n equations: 

a flog L) . = 0 

subject to the constraint that the bearing lines determined by 
the estimates must all pass through a common point. Because of 
this constraint, the number of independent estimates is reduced 
from n to 2. This is equivalent to stating that any two of 
the estimates determine a common point. The common point is the 
maximum likelihood estimate for the position of the object. 
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In order to impose the constraint on the estimates, the 
likelihood function will be reformulated in a rectangular 
coordinate system. To do this, first, consider the following 
definitions: 

u^ = r. ( 0 . - 0^) , = r^ (0^ - and w^ = r^ (0^ - 

where i = 1 to n is the observation number and the observations 
are from two or more stations. In these definitions, is the 

bearing of an initial estimate of the position of an object and 
r^ is its range from the station associated with the i ^ 
observation. It is also the radius of a circle that is centered 
on the station and passes through the initial estimate as shown 
in Figure 2. With , 0^ and in radians, u^ , v^ and Wj 

are the length of arcs on this circle. And, u^ = w^ - v- . In 
this relation, v. involves the unknown true bearing 0^ , but w- 
is known. 

In order to express v^ in a convenient form, consider a 
rectangular coordinate system whose origin is at the initial 
estimate, whose positive y-axis is north and whose positive 
x-axis is east. With x and y the coordinates of the object, 
to first order: v^ = x cos - y sin . This implies: 

0, = + (X cos - y sin / r^ 

And, to first order, this expresses the constraint on the 0^ 
and on their maximum likelihood estimates and it implies that the 
bearing lines determined by , 0^ and are parallel. 
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North 




Figure 2. In the program that implements the procedure, the 
origin of the xy-coordinate system is coincident with the common 
reference point rather than with the initial estimate and the 
coordinates of the initial estimate are x* and y*. The 
program is discussed later in the section. 
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In order to determine the estimates 0^ to first order, use 
the above relation to replace 0^ by the known values and 

the two unknown values x and y. Next, determine the maximum 
likelihood estimates x and y for the unknown coordinates x 
and y which, to first order, are the unknown coordinates of the 
object. That is, find the solutions to the equations determined 
by: 



3 flog L) 
5x 



= 0 



and 



x=x 

y=y 



3 flog L) 

ay 



= 0 



x=x 

y=y 



Evaluating the derivatives, the equations are: 



[ (w^ 


X cos )3i 


+ y 


sin 


^i) 


cos /3i] 


/ = 0 


[ (w^ 


X cos a, 


+ y 


sin 


^0 


sin )3i] 


/ a,. = 0 



where cr^ = r^ e^. These equations can also be written as: 

A X + B y = D 
B X + C y = E 



and their solutions as: 

X = (CD - BE) / (AC - B^) 

y = (AE - BD) / (AC - B^) 

where 

A = (cos^ ^i) / , B = (cos sin , 

C = SJ (sin^ ^-) / , D = (w^ cos 

E = 2j (w. sin /3.) / aj . 
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If the distribution of the random variables that determine x 
and y have a known distribution, then, in principle, a 
confidence region for an estimate can be constructed. With this 
in mind, the distribution of the estimators that determine x 
and y will be determined next. If X and Y represent the 
estimators, that is, the random variables that determine x and 
y / then, 

X = (Wi / o\) (B sin - C cos / (B^ - AC) 

Y = (W. / o]) (A sin - B cos / (B^ - AC) 

where W^ = r^ (6^ - ^,) . Since the 0^ are normally distributed 
random variables, the W^- are also. And, since the estimators X 
and Y are linear combinations of the , they have a 

bivariate normal distribution. 

If = 0, , for i = 1 to n, that is, if the initial 

estimate is the position of an object then, with E again the 
expected value operator, E (W J = 0 , since W^ = r^ (0^ - ^^) . In 

this case: 



— 2 


(1/a^ 


,) (B 


sin 




- c 


COS 




V(B^ - 


AC)^ 


ol = 2? 


(1/a^ 


i) (A 


sin 




- B 


cos 


^i) 


V(B^ - 


AC)^ 


^xy = 2? 




) (B 


sin 




- C 


cos 




(A sin 


- B COS - AC)^ 



Using the definitions for A , B , and C , that are given 
above, these equations can be written as; 
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a? = C / (AC - B^) 

a? = A / (AC - B^) 

aj- = B / (AC - B^) 

To the degree that the initial estimate is near the position of 

an object, these equations define approximations for the elements 
of the covariance matrix of the estimators X and Y. 

Since, in general, a-- will not be zero, the axes of the 
ellipse that bounds a minimum area confidence region for x and 
y will not be parallel to the xy-coordinate axes. To determine 
the orientation, translate the xy-coordinate system so that its 
origin is at the point determined by x and y . Then, rotate 
this coordinate system until its axes coincide with the axes of 
the bounding ellipse. Next, represent coordinates in this system 
by x' and y'. If the rotation angle is y , then the 
coordinates of a point in the two systems are related by the 
following general transformation equations: 

x' = X cos 7 - y sin 7 
y' = X sin 7 - y cos 7 

These transformation equations can be used to determine the 
elements of the covariance matrix of the distribution in the 
primed coordinate system from the covariance elements of the 
distribution matrix in an unprimed coordinate system, by using 
the following equations: 
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ol = 


4 


cos^ 


7 - 


2 a — 

xy 


sin 7 


cos 


7 + 


2 • 2 

oi sin^ 


oj. = 


4 


sin^ 


7 + 


2 a— 

xy 


sin 7 


cos 


7 + 


a? cos^ 


Q 

x> 

II 


(4 


- 


sin 


7 cos 


7 + 


a-- 

xy 


(cos^ 7 


- sin^ 



since X' and Y' are independent random variables when the axes 
of the coordinate system are parallel to the axes of the bounding 
ellipse, Oi'y' = 0 which implies that 7 satisfies the equation: 

tan 27 = / (a? - a?) . 

If the position of an object were known, and the estimated 
position were at the true position, then the probability p that 
an estimated position would be in a region bounded by an ellipse 
with axes 2ka*. and 2ka-. that is centered on the true position 
is given by p = 1 - exp(-kV2). This result can be obtained by 
integrating over the region bounded by the ellipse the bivariate 
normal density that determines X' and Y'. Note that k is a 
measure of the size of the region bounded by the ellipse and, for 
a probability of containment p , the corresponding value of k 

, ] /p 

IS determined by k = [-2 log(l - p) ) ] • These relations 

provide a basis for determining confidence regions for the 
estimates x and y of the position of an object that are given 
above . 

A program called PEST that is listed in Appendix 2 can be 
used to generate position estimates and corresponding confidence 
regions from bearing observations on or from two or more 
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stations. For a bearing on a station, 6 is the reciprocal of 
the observed bearing. In the program, the initial estimate is 
the point of intersection of the two bearing lines that are 
determined by the first two entered bearings with each from a 
different station. After the first two bearings have been 
entered, the program determines the coordinates x* and y* of 
the initial estimate by solving the following equations: 

X* sin ($2 - 0 ^) = [/Oj sin (a^ - ] sin $2 

- [ p 2 (®2 “ ^2)1 ^1 

y* sin ($2 - 61 ) = [yOj sin (ttj - 5 j) ] cos $2 

- [P2 sin (Q !2 - ^2)] ^1 

where ttj and Pj are the bearing and range of the first 
station and 02 and P 2 are the bearing and range of the sscond 
station from a common reference point. This method of finding 
the initial estimate is based on one that is used in Reference 3 . 
The point estimates d and p of the bearing and range of an 
object from a common reference point are shown in Figure 3. 

As an example, consider the data in Table 2. 



TABLE 2 



Station 

Number 


Observed 

Bearing 


Station 

Bearing 


Station 

Range 


Bearing 

Error 


1 


38° 


334° 


13500 


4° 


2 


324° 


50° 


11350 


3° 


3 


3° 


0° 


0 


4° 
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North 




East 



Figure 3. A confidence region bounded by an ellipse and a primed 
coordinate system in which is zero. The center of the 

ellipse and the origin of the primed and the unprimed coordinate 
systems are at the estimated position of an object. The 
direction of the major axis of the ellipse is 6 where 
0° < 5 < 180°. The rotation angle from the unprimed to the 
primed coordinate system is 7 . The estimates d and p are 
the bearing and range of the object from a common reference 
point . 
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The order of data entry in the program is determined by the 
station number. Note that the reference point is at the location 
of Station 3. If the distance is in meters, then, for the three 
observations, the estimates generated by the program for the 
bearing and range of the object from Station 3 are: a = 0.15° 

and p = 19554 meters. For k = 2, the lengths of the major and 
minor axes of the ellipse bounding the minimum area confidence 
region are: MJ = 3426 meters and MI = 2260 meters and the 

direction of its major axis is 12.48°. 

An extension of the above procedure is developed next. The 
extension accounts for cases where not all of the positions of 
the reference stations can be considered to be known. In these 
cases, if the station position were to become known, then, for a 
bearing line that is determined by an observation, the bearing 
line that is drawn from an assumed station position will be 
parallel to the bearing line that is drawn from the station 
position. The procedure is based on this relation and on the 
model with the addition of the following condition: The 

coordinates of an uncertain station position are determined by a 
distribution which is bivariate normal with a covariance matrix 
whose elements are known and a mean vector whose elements are the 
unknown coordinates of the station. This implies that for a 
bearing determined by an observation, the algebraic distance 
between the bearing line from an assumed station position and the 
bearing line from the station position is a random variable S 
that has a normal distribution with mean zero and a standard 
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deviation 0 ^ that is determined by the elements of the known 
covariance matrix. 

The geometry of the model is shown in Figure 4. There, 6 
is an observed bearing and the common origin of the coordinate 
systems is the assumed station position. For the unprimed 
system, the positive y-axis direction is north and the positive 
x-axis direction is east. For the primed system, the x'-axis 
is coincident with the major axes of an ellipse that bounds a 
confidence region for the assumed position. The direction of the 
major axis is 6 where 0° < 6 < 180°. For the double primed 
system, the y"-axis is coincident with the assumed bearing line. 

The effect of station position uncertainty is accounted for 
by the algebraic distance s between a bearing line from an 
assumed station position and the parallel bearing line from the 
station position. In effect, to first order, the following 
bearing lines are also parallel: the bearing line from a station 
to the initial estimate, the bearing line from the station to an 
object and a bearing line determined by an observation. Hence, 
the random variable = r^ (0, - 0^) is the algebraic distance 
between the observed bearing line and the bearing line from the 
station to an object. And Ug^ is a random 

variable that determines the algebraic distance between the 
observed bearing line from an uncertain station position and the 
bearing line from the station position to the object. Since U, 
and S are independent random variables, the standard deviation 
of Ug, is; ag, = [a? + . 
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assumed 
bearing line 




I 



X 



Figure 4. The geometry associated with the determination of the 
standard deviation . Here, 6 is an observed bearing and 
the common origin of the coordinate systems is the assumed 
station position. The positive y-axis direction is north and 
the positive x-axis direction is east. The x'-axis is 
coincident with the major axes of an ellipse that bounds a 
confidence region for the assumed station position, S is the 
direction of its major axis and the y"-axis is coincident v/ith 
the assumed location of the observed bearing line. 
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In the extended procedure, for a station of uncertain 
position, the standard deviation = r^e^ is replaced by: 






+ 



] 



1/2 



where 



= ol- sin^ (8-6) + CTy. cos^ (6 - 6) 

based on Figure 4. 

In Figure 4, note that ka^. and ka^. are the lengths of 
the major and minor axes of the bounding ellipse for some value 
of k. For an uncertain station position, the bearing and 

r^ the range of the initial estimate from the station position 
are determined by the assumed station position. These 
relationships are illustrated in Figure 5. 

In the development, although the position of one or more of 
the stations may be unknown, the position of a common reference 
point is known. This implies that the assumed station positions 
and the position of the initial estimate are also known. 
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object 



O 




Figure 5. In the development for the extended procedure, ^ 
r are the range and bearing of the initial estimate from the 
assumed station position. 



and 
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IV. A Position Estimate Based on Lines of Position 

The procedure that is described in this section is for 
determining a position estimate for an object from two or more 
lines of position that can be considered to be straight lines in 
the neighborhood of the object. The procedure is based on the 
following model: A line of position is a straight line. A line 

of position that is based on an observation is parallel to the 
true line of position. The algebraic distance u between a line 
of position that is based on an observation and the true line of 
position is the value of a normally distributed random variable 
with zero mean and known standard deviation. 

Labeling lines of position by the index i , the model 
implies that the algebraic distance Wj between the i th line of 
position based on an observation and a parallel line through an 
initial estimate of an unknown position is related to the 
algebraic distance u^ between it and the true line of position 
by: u^ = w^- v^ where v^ is the algebraic distance between the 

i th true line of position and the parallel line through the 
initial estimate. Now, consider a rectangular coordinate system 
whose origin is at the initial estimate, whose positive y-axis 
is north and whose positive x-axis is east. With x and y 
the coordinates of the unknown position, to first order: 

v^ = X cos - y sin 

where is the direction of the i th true line of position and 

0° < < 180°. By referring to the development in Section III, 
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note that the estimates x and y that are developed there also 
apply here. However, in this development the standard deviation 
a of the random variable that determines the distance u must be 
known for each of the observed bearing lines. 

Figure 6 illustrates the geometry associated with the 
development. 
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Figure 6 . The geometry for three lines of position. The 
direction of the lines are , 02 ^ 3 * With the initial 

estimate of the unknown position at the intersection of line 1 
and line 2 , both Wj and W 2 are equal to zero. This implies 
that the final estimate of the position and the initial estimate 
of the position are equal. If the three lines of position v.’ere 
obtained by sextant observations, the assumed position should be 
the reference position. In this case, the angles Qj , and 

<13 would be the azimuth angles. 
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Appendix 1. The Program COMP 



COMP is a program with which to implement the procedure 



developed in Section II for generating a composite estimate of 



the position of an object from two or more independent estimates 



of the position of the object. The independent estimates must be 



specified in terms of coordinates in a rectangular coordinate 



system whose positive y-axis is north and positive x-axis is 



east. In addition, a confidence region for each estimate must be 



specified in terms of the direction 6 of the major axis of a 



bounding ellipse and either its size k or its containment 



probability p (confidence) . 



The program listing follows: 



10 CLS : PI = 4 * ATN(l) 

20 INPUT "print data yes/no (y/n)”; A$ 

30 IF A$ = "N" OR A$ = "n" THEN P$ = "N"; GOTO 60 
40 IF A$ = "Y" OR A$ = "y" THEN P$ = "Y" ELSE GOTO 20 

50 LPRINT : LPRINT : LPRINT : LPRINT "Composite Position Estimation Program: COMP. BAS" 

60 PRINT : INPUT "estimates to be combined": NO 
70 IF P$ = "N" THEN GOTO 90 

80 LPRINT : LPRINT "estimates to be combined ", NO 

90 INPUT "elliptical confidence regions by size or probability (s/p)"; F$ 

100 IF F$ = "S" OR F$ = "s" OR F$ = "P" OR F$ = "p" THEN GOTO 110 ELSE GOTO 90 

no IF F$ = "P" OR F$ = "p" THEN B$ = "probability" ELSE B$ = "size" 

120 FOR 10 = 1 TO NO 

130 PRINT : INPUT "position estimate x-coordinate"; X 

140 INPUT "position estimate y-coordinate"; Y 

150 INPUT "containment ellipse major axis direction"; 00 

160 IF 00 >= 180 OR 00 < 0 THEN GOTO 150 

170 INPUT "ellipse major axis length"; MA: SMA = MA / 2 

180 INPUT "ellipse minor axis length": MI: SMI = MI / 2 

190 IF F$ = "P" OR F$ = "p" THEN GOTO 220 

200 INPUT "containment ellipse size"; SO: PO = 1 - EXP(-S0 * SO / 2) 

210 PRINT "containment probability = "; PO: GOTO 240 

220 INPUT "containment probability"; PO: SO = SQR(-2 * L0G(1 - PO)) 

230 PRINT "containment ellipse size = "; SO 
240 IF P$ = "N" THEN GOTO 320 

250 LPRINT : LPRINT : LPRINT "position estimate x-coordinate", X 

260 LPRINT "position estimate y-coordinate", Y 

270 LPRINT "ellipse major axis direction", 00 

280 LPRINT "ellipse major axis length ", MA 

290 LPRINT "ellipse minor axis length ", MI 

300 LPRINT "ellipse size ", SO 

310 LPRINT "containment probability ", PO 

320 SX = SMA / SO: SY = SMI / SO 

330 00 = 00 * PI / 180 

340 U = SX * SX * SIN(DO) * SIN(OO) + SY * SY * COS(OO) * COS(OO) 

350 V = SX * SX * COS(OO) * COS(OO) + SY * SY * SIN(OO) * SIN(OO) 

360 W = (SX * SX - SY * SY) * SIN(OO) * COS(OO) 

370Z=U*V-W*W:A=A+V/Z:B=B-W/Z:C=C+U/Z:0=0+V/Z*X-W/Z*Y 
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380E=E-W/Z*X+U/Z*Y:F=F+V/Z*V/Z*U:G=G+W/Z*W/Z* 

390 H = H- V/ Z*W/Z*U: I = I+ W/ Z*W/Z*V:J = J + U/ Z*U/Z’^V 

400K=K-W/Z*U/Z*V:L=L-V/Z*W/Z*W 

410M = M+(V/Z’"U/Z + W/ Z*W/Z)*W:N = N- W/ Z*U/Z*W 

420 NEXT 10 

430Q=A*C-B*B:X=(C*D-B*E)/Q: Y=(A^E-B*D)/Q 

440 R=(C*C*(F+I + 2*L)-2*C*B*(H + K + M)+B*B^(G + J + 2’^N)) 

450 S=(A*A*(G + J + 2*N)-2’^A*B*(H + K + M) + B*B*(F-fI+2*D) 

460 T = ((A * C + B * B) * (H + K + M) - C ^ B * (F -K I + 2 * L) - B A ^ (G + J + 2 

470 CLS : PRINT ’’composite estimate x-coordinate = X 

480 PRINT ’’composite estimate y-coordinate = Y 

490 IF P$ = ”N” THEN GOTO 520 

500 LPRINT : LPRINT : LPRINT : LPRINT ’’composite estimate x-coordinate”, X 
510 LPRINT ’’composite estimate y-coordinate”, Y 
520 IF S = R THEN A = 0: GOTO 600: REM a circle 

530 A = .5 * ATN(2 * T / (S - R)) 

540 R1 = R * COS(A) * COS(A) - 2 * T * COS(A) * SIN(A) + S * SIN(A) * SIN(A) 

550 SI = R * SIN(A) * SIN(A) + 2 * T * COS(A) * SIN(A) + S * COS(A) * COS(A) 

560 A = A * 180 / PI 

570 IF R1 > SI THEN A = A + 90: GOTO 590 



580 RT = R1 
590 IF A < 
600 PRINT 
610 IF G$ ^ 
620 IF G$ ^ 
630 PRINT 
640 PRINT 
650 PRINT 
660 PRINT 
670 PRINT 
680 PRINT 
690 PRINT 
700 IF P$ 
710 LPRINT 
720 LPRINT 
730 LPRINT 
740 LPRINT 
750 LPRINT 
760 PRINT 
770 IF 6$ = 
780 IF G$ = 
790 END 



R1 = SI: SI = RT 
0 THEN A = A + 180 

PRINT : INPUT ’’confidence region by size or probability (s/p)”; G$ 

”S” OR G$ = ”s” THEN GOTO 630 

”P” OR G$ = ”p” THEN GOTO 650 ELSE GOTO 600 

INPUT ’’containment ellipse size”; SO: PO = 1 - EXP(-S0 * SO / 2) 
PRINT "containment probability = ”; PO: GOTO 670 



SO = SQR(-2 
SO 

* SQR(Rl) 



2 * SO * SQR(Sl) 



SO 



INPUT "containment probability"; PO: 

PRINT "containment ellipse size = ”; 
ellipse major axis length = ”; 2 * SO 
ellipse major axis direction = ’’; A 
ellipse minor axis length = 

: "N" THEN GOTO 760 
: LPRINT ’’containment ellipse size 

"containment probability ”, 

"ellipse major axis length ”, 

"ellipse major axis direction ”, 

"ellipse minor axis length 
PRINT : INPUT "continue or quit (c/q) 

”C" OR G$ = "c” THEN GOTO 600 

”Q” OR G$ = ”q” THEN GOTO 790 ELSE GOTO 770 



L0G(1 - PO)) 



PO 
2 * 
A 

2 * 



SO 



SO 



G$ 



SQR(Rl) 

SQR(Sl) 



(Q " Q) 

(Q * 0) 

N)) / (Q * Q) 
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Appendix 2. The Program PEST 



PEST is a program with which to implement the procedure 



developed in Section III for determining a position estimate for 



an object from a set of observed bearings from or on two or more 



reference points of known positions. In the development, the 



reference points are referred to as stations. The bearing error 



associated with an observation must be specified for each 



station. In addition, the position of each station must be 



specified in terms of its bearing and range from a single 



reference point which may be one of the reference points. The 



program generates an estimate for the position of an object and a 



confidence region for the estimate which is specified in terms of 
the direction 5 of the major axis of a bounding ellipse and 



either its size k or its containment probability p 



(confidence) . 



The program listing follows: 



10 as : PI = 4 * ATN(l): DIM A(7); I = 0 
20 INPUT "print data (y/n)"; A$ 

30 IF A$ = "N" OR A$ = "n" THEN P$ = "N": GOTO 60 

40 IF A$ = "Y" OR A$ = "y" THEN P$ = ”Y" ELSE GOTO 20 

50 LPRINT: LPRINT: LPRINT : LPRINT "Position Estimation Program: PEST. BAS" 
60 INPUT "bearings on object or bearings on stations (o/s)": B$ 

70 IF B$ = "0" OR B$ = "o" THEN K = 0: GOTO 120 

80 IF B$ = "S" OR B$ = "s" THEN K = 1 ELSE GOTO 60 

90 IF P$ = "N" THEN GOTO 120 

100 IF K = 0 THEN B$ = " object" ELSE B$ = " stations" 
no LPRINT : LPRINT "bearings on" + B$: LPRINT 
120 CLS : PRINT "input number " + STR$(I + 1) 

130 INPUT "observed bearing in decimal degrees": P 
140 INPUT "station bearing in decimal degrees": Q 
150 INPUT "station range": R 

160 INPUT "bearing error in decimal degrees": BE 
170 IF P$ = "N" THEN GOTO 230 

180 LPRINT : LPRINT "input number " + STR$(I + 1) 

190 LPRINT "observed bearing in decimal degrees". P 
200 LPRINT "station bearing in decimal degrees", Q 
210 LPRINT "station range in chosen units", R 
220 LPRINT "bearing error in decimal degrees", BE 
230 IF K = 1 THEN P = P + 180 

240 P = P * PI / 180: Q = Q * PI / 180: BE = BE * PI / 180 
250 IF I = 2 THEN GOTO 350 

260 I = I + 1: A(I - 1) = P: A(I + 1) = Q: A(I + 3) = R: A(I + 5) = BE 
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270 IF I = 1 THEN GOTO 120 

280 X = A(4) * SIN(A(2) - A(0)): Y = A(5) * SIN(A(3) - A(l)); Z = SIN(A(1) - A(0)) 

290 IF Z = 0 THEN GOTO 900 

300 U = (X * SIN(A(D) - Y * SIN(A(0))) / Z; V = (X * C0S(A(1)) - Y * C0S(A(0))) / Z 
310 FOR J = 0 TO 1 

320 P = A(J): Q = A(J + 2): R = A(J + 4): BE = A(J + 6): GOSUB 940 
330 NEXT J 
340 GOTO 360 
350 GOSUB 940 

360 CIS : INPUT "continue data input or generate an estimate (c/g)"; C$ 

370 IF C$ = "G" OR C$ = "g" THEN GOTO 390 

380 IF C$ = "C" OR C$ = "c" THEN GOTO 120 ELSE GOTO 360 

390 CLS 

400 F = (B * B - A * C): IF F = 0 THEN GOTO 900 

410 XI = U + (B * E - C D) / F: Y1 = V + (A * E - B * D) / F: GOSUB 1000 

420 R2 = Rl: B2 = B1 * 180 / PI : IF K = 1 THEN B2 = B2 - 180: IF B2 < 0 THEN B2 = 360 B2 

430 T = SGN(B) * PI / 4: IF A = C THEN GOTO 450 

440 T = .5 * ATN(2 * B / (A - C)) 

450 G = (C * COS(T) * COS(T) - 2 * B * COS(T) * SIN(T) + A * SIN(T) * SIN(T)) / (-F): G = SQR(G) 

460 H = (C * SIN(T) SIN(T) + 2 * B * COS(T) * SIN(T) + A * COS(T) * COS(T)) / (-F): H = SQR(H) 

470 IF H >= G GOTO 490 

480 Z = H: H = G: G = Z: T = T + PI / 2 

490 CLS 

500 PRINT "bearing estimate = B2 
510 PRINT "range estimate = "; R2: PRINT 
520 IF P$ = "N" THEN GOTO 550 

530 LPRINT : LPRINT : LPRINT "bearing estimate", B2 
540 LPRINT "range estimate", R2 

550 INPUT "confidence region by size or probability (s/p)"; C$ 

560 IF C$ = "S" OR C$ - "s" THEN GOTO 580 

570 IF C$ = "P" OR C$ = "p" THEN GOTO 650 ELSE GOTO 550 

580 INPUT "confidence region size"; S 

590 IF S <= 0 THEN GOTO 580 

600 CP = 1 “ EXP(-S S / 2) 

610 PRINT : PRINT "containment probability = "; CP 

620 IF P$ = "N" THEN GOTO 720 

630 LPRINT : LPRINT "confidence region size", S 

640 LPRINT "containment probability", CP: GOTO 720 

650 INPUT "containment probability (confidence)"; CP 

660 IF CP >= 1 OR CP <= 0 THEN GOTO 650 

670 S = SQR(-2 * L0G(1 - CP)) 

680 PRINT : PRINT "confidence region size = S 

690 IF P$ = "N" THEN GOTO 720 

700 LPRINT : LPRINT "containment probability", CP 

710 LPRINT "confidence region size", S 

720 X = 2 * S * H 

730 PRINT "major axis length = X 

740 N = T * 180 / PI: IF N < 0 THEN N = N + 180 

750 PRINT "major axis direction = N 

760 Y = 2 * S * G 

770 PRINT "minor axis length = "; Y 

780 PRINT "containment ellipse area = "; PI * X * Y / 4 

790 IF P$ = "N" THEN GOTO 840 

800 LPRINT "major axis length", X 

810 LPRINT "major axis direction", N 

820 LPRINT "minor axis length", Y 

830 LPRINT "containment ellipse area", PI * X * Y / 4 
840 PRINT : PRINT ; INPUT "continue or quit (c/q)"; C$ 

850 IF C$ = "C" OR C$ = "c" THEN GOTO 870 

860 IF C$ = "Q" OR C$ = "q" THEN END ELSE GOTO 840 

870 CLS : INPUT "continue data input or generate an estimate (c/g)"; C$ 

880 IF C$ = "G" OR C$ = "g" THEN PRINT : GOTO 500 

890 IF C$ = "C" OR C$ = "c" THEN GOTO 120 ELSE GOTO 430 

900 PRINT "no solution" 

910 IF P$ = "N" THEN GOTO 930 
920 LPRINT : LPRINT "no solution" 

930 END 

940 XI = U - R * SIN(Q): Y1 = V - R * COS(Q): GOSUB 1000 
950 W = P - Bl: L = Rl * BE: IF L = 0 THEN GOTO 900 
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960 6 = COS(Bl) / L: H = SIN(Bl) / L 

970 IF U >= PI THEN W = W - 2 * PI: GOTO 990 

980 IF U <= -PI THEN W = W + 2 * PI 

990W=W/BE:A=G*G+A:B=G*H+B:C=H*H+C:D=W*G+D:E=U*H+E: RETURN 
1000 R1 = SQR(X1 * XI + Y1 * Yl): IF R1 = 0 THEN B1 = 0: RETURN 

1010 IF ABS(X1 / Rl) = 1 THEN Ml = SGN(Xl) * PI / 2 ELSE Ml = ATN(X1 / R1 / SQR(1 - XI * XI / R1 / Rl)) 

1020 IF ABS(Y1 / Rl) = 1 THEN B1 = (1 - SGN(Yl)) * PI / 2 ELSE B1 = PI / 2 - ATN(Y1 / Rl / SQR(1 - Yl * Yl / 
Rl / Rl)) 

1030 IF Ml < 0 THEN B1 = 2 * PI - B1 
1040 RETURN 
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